NASA Technical Paper 1910 


NASA 
TP 
1910 
c. 1 


V iscous-Inviscid Interaction 
Method Including Wake 
Effects for Three-Dimensional 
Wing-Body Configurations 


.mTECHm»n= L 
K1RTI.AND AFE a 2 

cr pr 

Ln -n 

sss m 










NASA Technical Paper 1910 


TECH LIBRARY KAFB, NM 



□□b7bfi5 


Viscous-Inviscid Interaction 
Method Including Wake 
Effects for Three-Dimensional 
Wing-Body Configurations 


Craig L. Streett 
Langley Research Center 
Hampton, Virginia 


IWNSA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1981 




SUMMARY 


An existing three-dimensional compressible integral boundary-layer method 
was modified to account for mean dilatation effects, to model transition prop- 
erly, and to provide better numerical stability near computational boundaries- 
Results of this method were compared with those from a three-dimensional finite- 
difference boundary-layer method on a difficult test case. An interaction pro- 
cedure was developed to couple this integral method with a number of wing-alone 
and wing-body transonic potential codes to account for viscous effects- A 
locally two-dimensional wake model, including thickness and curvature effects, 
was developed and incorporated into this interaction procedure. Results from 
this procedure were compared with experimental data and results from previous 
procedures for test cases where viscous effects were large- 


INTRODUCTIC»I 

It is now generally recognized that any analysis procedure for canputing 
transonic flew over wing and wing-body configurations must include some form of 
viscous correction on the wing to give reasonable agreement with experiment. 

A number of researchers have shown good agreement without viscous correction 
through the use of lift-matching; however, such techniques require a priori 
knowledge of the overall lift, rather than it being a computational result. 
Ideally, one would like to compute flows using only known geometric quantities 
and flow conditions. Some form of viscous correction thus appears to be 
required. 

Techniques for solving the Navier -Stokes equations, either in full form or 
in some reduced form (parabolized, velocity split, etc.), are still in their 
infancy. In addition, Navier -Stokes computations for a transonic wing at 
Reynolds numbers of interest are still beyond the capability of present comput- 
ers. Alternatively, potential- flow and boundary-layer computational techniques 
are mature, and the viscous-inviscid interaction philosophy of computing flow 
fields with viscous effects has been applied successfully in the past to prob- 
lems of interest. 


Viscous-inviscid interaction consists of iteratively solving the boundary- 
layer and potential (or other inviscid flow) equations, coupled through the con- 
cept of boundary- layer displacement thickness. If the Navier-Stokes equations 
are written in terms of matched asymptotic expansions near and far from a wall, 
it can be shown that to 0(Nj^e“^/^) the pressure on the wall is determined from 
the solution of the outer (inviscid) equations, with the flow tangency boundary 
condition displaced from the wall by a distance determined from the inner 
( boundary- layer ) equations. This viscous-inviscid interaction technique produces 
a uniformly valid solution, with the exception of strong interaction regions, for 
example, near the trailing edge, a normal shock, or the wing tip. 



Illllllllllllllllllllllllll 


Melnik, Chow, and Mead (ref. 1) have shown in two dimensions that this 
interaction of the first-order equations involves a singularity in the pres- 
sure gradient and streamline curvature at the trailing edge of an airfoil. It 
has been demonstrated computationally, however, that iterating the boundary layer 
and potential equations until convergence produces a self-consistent solution in 
which this singularity does not appear. It is shown in reference 2 that although 
this self-consistent solution is not singular, it is inconsistent with the normal 
manentum equation in a small region near the trailing edge. Correction to the 
inviscid boundary conditions near the trailing edge in two dimensions is devel- 
oped in reference 1 to account for this inconsistency; a similar development for 
three-dimensional wings with spanwise flow at the trailing edge is not available. 
Although self-consistent interaction solutions are not rigorously accurate to 
first order near the trailing edge, such solutions do show good agreement with 
experiment. Such a self-consistent development is described in this paper. 

There are several boundary-layer schemes which may be applied to an inter- 
action procedure. Several researchers (refs. 3 to 5) have demonstrated fair 
agreement with experiment by using two-dimensional boundary-layer calculations 
applied in chordwise strips over the wing span. However, experiments with 
modern transport wings with moderate sweep, utilizing supercritical airfoil 
sections, show large amounts of spanwise flow in the lower surface cove region. 
This flow tends to relieve the mass buildup in the lower surface boundary layer 
near the root and increase it near the tip. Strip theory would then predict 
more viscous decambering near the root and less near the tip than a method which 
accounts for this cross flow. When this effect is significant, a fully three- 
dimensional viscous correction should provide improved accuracy compared with 
two-dimensional str i p theory . 

A number of finite-difference procedures have been developed to compute wing 
boundary layers in three dimensions, such as the one described in reference 6. 
These methods still require considerable empiricism in the form of a turbulence 
model, and they are expensive to run. These codes also tend to be quite diffi- 
cult to use in an interactive scheme, as primary flow separation and inflow at a 
boundary are often too destabilizing ccxnputationally for the procedure to with- 
stand. Thus the boundary-layer procedure may fail to give a solution in sane 
region of the wing so that further global iterations are impossible. This situ- 
ation can occur even when the predicted boundary layer at convergence is 
well-behaved. 

Three-dimensional integral methods are less complex and more reliable 
than the finite-difference methods but present a more physically reasonable 
description of the three-dimensional boundary layer than the two-dimensional 
strip methods. In general, integral methods are characterized by replacing 
detail velocity profile calculations with assumptions of general profile shapes 
or families. These profile families are characterized by one or more integral 
parameters, such as manentum thickness and shape factor. In the three- 
dimensional case, the boundary layer is split into streamwise (i.e., parallel 
to the external streamline) and cross-flow (perpendicular to the external 
streamline) components, as was suggested by Coles (ref. 7). These profile 
assumptions are substituted into the governing equations, which are in turn 
integrated across the boundary layer in the wing surface normal direction. 

Thus, the three-dimensional boundary- layer equations are converted into a 
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set of two-dimensional equations in the planform plane for the unknown integral 
parameters. 

The three-dimensional integral boundary-layer method used in this work 
was based on a code developed by the German aerospace firm, Dornier G.m.b.H. 

(ref. 8) • A description of the method, along with the modifications incorpo- 
rated into it during this study, appears in a later section. Results from this 
code are compared with those of a three-dimensional finite-difference boundary- 
layer code. The interaction of this code with two transonic full- potential wing- 
alone codes and one wing-body code is described, and results for two transport 
wing-body test cases with supercritical airfoil sections are shown. For one 
test case, overall viscous effects, body effects, and influence of different 
wake treatments are described. 

It is shown in this paper that, 1:^ including the effects of wing surface 
boundary layer and both wake displacement and curvature, wing lift distribution 
and shock position are predicted well. However, since the strong viscous- 
inviscid interaction of the trailing-edge region is not modeled properly, 
trailing-edge pressures are generally predicted too high. The drag predicted 
by this method, which accounts only for weak interaction, thus is in error. 

It is shown in reference 1 for two-dimensional airfoil flow that the trailing- 
edge strong interaction region must be accounted for to obtain accurate drag 
prediction. 


SYMBOLS 

ag speed of sound at edge of boundary layer 

ao stagnation speed of sound 

A aspect ratio 

b semispan 

c local chord 

Cl local section lift coefficient 

Cp pressure coefficient 

DFVLR Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt E.V. 
G reduced potential used in inviscid codes 

H streamwise shape factor, = Oi/0-| 

M Mach number 

^Re Reynolds number per unit length 

^Re,c Reynolds number based on chord 
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R normalized lift reduction due to viscous effects 

Sy arc length along wake 

t maximum section thickness 

u local boundary-layer velocity in external streamwise direction 

U external velocity magnitude 

V local boundary-layer velocity normal to external streamwise direction 

X chor dwi se coordinate 

y spanwise coordinate 

z j^ysical normal coordinate 

Z transformed normal coordinate 

3 angle between external streamline and wall shear direction 

wake vorticity 

6* total boundary-layer displacement thickness 
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streamwise displacement thickness. 
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cross- flow displacement thickness. 



dz 
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* 

w 


wake displacement thickness. 




dz 


A transformed boundary-layer thickness 

e local flow angle in constant-span plane 

n normalized spanwise coordinate, = y/b 
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streamwise momentum thickness, = 
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wake momentum thickness, = I — - (1 « - 
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k:„ wake curvature 

^c/4 quarter-chord sweep angle 

P local density 

Pg density at edge of boundary layer or wake 

Pq stagnation density 


BOUNDARY- LAYER METHOD AND VERIFICATION 

The boundary-layer method used in this study was obtained from Dornier 
G.m.b.H.; it is based on the laminar scheme developed in reference 8 and the 
turbulent method described in reference 9, with a number of modifications. 
Details of the development of these schemes is not given herein; rather, a 
general description and discussion of assumptions and their effects follow. 


Boundary-Layer Analysis 

For purposes of analysis, the boundary layer is decomposed into components 
parallel and perpendicular to the external streamline. A basic schematic of 
this decomposition of the boundary layer is shown in figure 1 . The situation 
shown is that most commonly encountered, in which the cross-flow velocity does 
not change sign across the boundary layer. 

The laminar analysis used in this study is basically unchanged from that 
developed in reference 8. The streamwise velocity profile is assumed to be of 
the Faulkner-Skan family of similarity profiles. The cross-flow profile is gen- 
erated from a linear combination of these similarity profiles; this combination 
is capable of generating profiles with cross-flow velocity crossover. These 
profiles are defined for incompressible flow; the Stewartson transformation is 
used to scale the normal coordinate: 

Pa^ 

dZ = dz 

^ 0^0 
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By this transformation, the compressible velocity profile may be represented by 
the incompressible form, as a function of the transformed variable Z. 

These profile assumptions are substituted into the two boundary-layer 
momentum equations and their corresponding moment of momentum relations. The 
moment relations are often used in integral boundary-layer methods (ref. 10); 
they are independent relations formed by multiplying the momentum equations by 
the distance from the wall. These four equations are then integrated from the 
wall to the boundary-layer edge. The resultant equations form a system of four 
coupled partial differential equations in two dimensions on the surface of the 
wing for four unknowns; the four unknowns are parameters related to the stream- 
wise momentum thickness, the wall shear, the cross- flow displacement thickness, 
and the boundary- layer thickness. The streamwise momentum thickness is defined 
as 


oo 



and the cross-flow displacement thickness as 




V 

- dz 
U 


In the turbulent case, a simple power-law profile is assumed for the 
streamwise direction in the form 


u 

U 



n(H) 


where Z is the transformed normal coordinate, and the exponent n is a func- 
tion of the streamwise shape factor H. The cross-flow profile used was sug- 
gested by Mager (ref. 11) as follows: 


V 

U 




tan 6 
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where 3 is the angle between the external streamline and the wall shear direc- 
tion. Note that this relation cannot predict a cross- flow profile with velocity 
crossover . 

These profiles are substituted into the boundary-layer momentum equations 
euid integrated, as before; a modified Ludwieg-Tillmann relation is used to eval- 
uate the wall shear (ref. 12). To close the system of equations, an integral 
continuity equation is formed, in which the change in mass flew in the boundary 
layer is explicitly identified. This change in mass flow, the so-called entrain- 
ment coefficient, is related to the other integral quantities through the lag- 
entrainment concept of Green et al. (ref. 1 3) . Thus, in the turbulent case, a 
set of three coupled first-order PDE's is generated for the streamwise momentum 
thickness, the shape factor, and the wall shear angle. The use of the lag- 
entrainment relations, rather than a local entrainment equation is the difference 
between the method used in this study and that described in reference 9. 

For a three-dimensional boundary layer, the displacement thickness must be 
computed as a solution to a PDE on the surface, rather than being computed from 
a local integral as in the two-dimensional case. In both the laminar and turbu- 
lent methods, therefore, an a posteriori solution of the integrated continuity 
equation yields the displacement thickness. 

Rather than using a streamline coordinate system, as is usually done in 
three-dimensional boundary-layer calculations, a simple wing constant-chord 
constant-span system is used. Although nonorthogonal, this system alleviates 
the difficulty of having to generate a complex new coordinate system for each 
viscous-inviscid Iteration. The additional terms from the nonorthogonal trans- 
formation are not computationally time-consuming because the systems of equations 
are first order; therefore, this coordinate system is considered more efficient 
than streamline coordinates for the integral methods. A schematic representa- 
tion of the coordinate system used is given in figure 2. 

The three-dimensional boundary-layer equations display hyperbolic proper- 
ties in the planform plane; that is, at any point on the plane, the equations 
have a dexnain of dependence which must be honored in any numerical solution 
scheme. The integral equations have the same property. It is shown in refer- 
ence 14 that the characteristics of the turbulent integral equations always lie 
between the external streamline and the wall shear line, giving a convenient, 
conservative approximation to the exact domain of dependence. 

An explicit finite-difference scheme is used to march the solution of the 
coupled system of equations in the chordwise direction. Differencing of (quanti- 
ties in the y-direction (spanwise) must be upwind to properly account for the 
(kimain of (3ependence. That is, if both the external streamline and wall shear 
line extend to the same side of the chord line frexn a computational point, the 
spanwise difference must be taken as one-sided from the opposite side. If the 
external streamline and wall shear line lie on opposite sides of the chord line, 
a central difference is used for the spanwise derivative. In addition, the 
step size in the x-direction (chordwise) is limited by the Courant-Frledricks- 
Lewy (CFL) condition, a necessary condition for the stability of an approximate 
solution scheme for a hyperbolic equation (ref. 15). 
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The spanwise differencing scheme of reference 14 must be modified at the 
wing root and tipr when flow into the computational region occurs at these 
boundaries, in order to make the procedure more reliable. Strictly, initial 
conditions must be provided at these locations. None being available, condi- 
tions along the boundaries must be assumed in these cases if the numerical 
scheme is to be stable. At the wing root, a plane of symmetry is assumed. Thus, 
on this chord line, the cross-flcw velocity is set to zero, as are all spanwise 
derivatives. Such an assumption is at least a reasonable model of the physical 
situation. During the course of the present study, this wing-root boundary 
condition was also approximated using infinite swept-wing assumptions, that 
is, spanwise derivatives only set to zero. This boundary condition, however, 
was not found to be strong enough to be reliable for general use. 

The wing-tip boundary is considerably more difficult. For this study, the 
assumption of zero spanwise derivatives on the tip chord line is used when flow 
into the region along this boundary is signaled. This is equivalent to assuming 
locally that the tip is part of an infinite swept wing. This assumption was 
found to be the weakest condition of those examined that was always stable and 
gave the least disturbance across the span. The inviscid method does not 
accurately model the outer flow in the tip region because of assumptions made 
concerning the shape of the trailing vortex sheet; therefore, the boundary-layer 
assumption made in the tip region is not considered restrictive. 

The initial condition along the leading edge is provided by an approximate 
attachment line analysis, such as that contained in reference 16. The attach- 
ment line is assumed to be along the wing leading edge. Transition location 
is input, fixed along a span line at that fraction chord. The transition 
model was modified frcm the original boundary-layer code obtained from Dornier 
to enforce continuity of the displacement thickness at transition. There is no 
inconsistency in this assumption as there would be in many two-dimensional inte- 
gral methods because 6* is not an integral parameter of either the laminar or 
the turbulent method but is developed from a separate integration of the con- 
tinuity equation. 

Separation is signaled in the streamwise direction when integral parameters 
reach critical values; in the laminar case, when the transformed wall shear 
reaches zero; and in the turbulent case, when the streamwise shape factor 
reaches 2.4. No cross- flow separation model is used. When laminar separation 
is detected along a chord line, transition is assumed over the entire span. If 
separation is signaled at a point in a turbulent region, the streamwise shape 
factor, skin friction, and wall shear angle are frozen thereafter along that 
chord line, and spanwise derivatives are set to zero. A one- dimensional inte- 
gration of the integral equations is thereafter performed on this chord line. 

This is considered an extrapolation technique; it should be understood that the 
integral equations with the given profile assumptions are invalid in large 
regions of separation. However, this method is assumed to give at least a rea- 
sonable solution in this situation and is required since such separated regions 
often appear only during the early iterations of the viscous- inviscid interac- 
tion. If the boundary-layer method were to stop at separation, the iteration 
procedure could not continue. If the interaction scheme converges to a solution 
with large separation regions remaining, the underlying assumption of the exist- 
ence of a thin viscous region is invalid. 
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Modifications were also made to the lag- entrainment calculations in order 
to improve the predictions. First, the streamwise pressure gradient was 

correctly incorporated into the equations which were originally developed for 
the two-dimensional case. Second, a modification to the dissipation length 
scale, due to dilatation of the compressible flow, was applied (ref. 13). The 
effects of these modifications are seen only in regions of severe pressure gradi- 
ent or Icurge cross flow. 


Verification of Boundary-Layer Method 

In order to gain confidence in the integral boundary-layer method, a rather 
severe test case was run for comparison with results from the three-dimensional 
finite- difference boundary-layer method described in reference 6. The velocity 
data input to both boundcury-layer codes was a final converged interacted solu- 
tion for the flow about a swept, tapered transport wing with significant aft 
camber. The interaction procedure used to produce this solution is the topic 
of later sections of this report; details of this procedure are not important 
here. A converged solution was used as input since the finite-difference code 
of reference 6 was unable to produce a boundary-layer solution over so much of 
the wing surface that comparisons were meaningless when the inviscid, uninter- 
acted velocity distribution was input. 

Figure 3 shows the difference in the displacement thickness predictions 
of the original Dornier three-dimensional integral boundary-layer code and the 
code modified for use in this study. Differences are seen near transition 
(which for this case occurred at 0.17 c), after the shock near midchord on the 
upper surface, at the upper surface trailing edge, and in the cove region on 
the lower surface. 

A comparison of the displacement thickness predictions of the modified 
integral method and of the finite-difference scheme of reference 6 is given in 
figure 4. The span station at which this comparison is made is near midspan 
since assumptions made at the wing root and tip are different for the two meth- 
ods. The agreement shown at this midspan station is typical for most of the 
wing. The largest percent difference is just aft of the midchord shock. The 
two methods also agree well on the location of the small upper surface trailing- 
edge separation region and on the wall shear angle in the lower surface cove 
region. Even the skin-friction coefficient predicted by the integral method is 
within 10 percent of that predicted by the finite-difference code, despite the 
very crude skin-friction model used in the integral method. The finite- 
difference code required almost 3200 CP seconds on the Control Data CYBER 1 73 
Computer System to produce a boundary-layer solution for this test case, whereas 
the integral code required less than one-thirtieth of this time and used about 
one-half the core storage. 

Results of the three-dimensional integral boundary-layer method were also 
compared with experimental data. Although three-dimensional test cases were 
desired, cases with adequate coverage of velocity data could not be found. In 
principle, the two momentum equations and the continuity equation could be 
solved in the planform plane at the edge of the boundary layer to derive the 
external velocity distribution given the surface pressure. However, such a 
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scheme requires fixing a set of rather arbitrary boundary conditions, especially 
at the critical tip region which can be an inflow boundary for the upper surface. 
In addition, the equations are first order and were found to be quite sensitive 
to the applied boundary conditions; this scheme was thus dropped after some 
effort, and some two-dimensional test cases were run. 

Perhaps the most representative and difficult two-dimensional test case 
was that of the upper surface boundary layer on the CAST 7 airfoil as tested 
in the DFVLR 1 x T meter tunnel (ref. 17). The case included a moderately 
strong shock and some separation at the trailing edge. The test and calculation 
were run at a Mach number of 0.76 and a Reynolds number of 2.4 x 10^. Calcula- 
tion was begun from measured momentum thickness and shape factor at 30 percent 
chord. Figure 5 shows the measured pressure distribution and the measured and 
computed momentum thickness distributions. The rapid thickening of the boundary 
layer in the high adverse pressure gradient region of the shock is predicted 
well, as is the thickness at the trailing edge. Proper prediction at the trail- 
ing edge is quite important, as lift is greatly influenced by the trailing-edge 
angle of the displacement body. Note the slight postshock disagreement due, 
perhaps, to the lack of a proper shock — boundary-layer interaction model. How- 
ever, since the experiment shows the thickness reducing in a region of zero 
pressure gradient, one must also suspect some three-dimensional relief in the 
experiment, perhaps the shock curving back near the sidewalls. 

Although the physical modeling of the boundary-layer flow is far more 
accurate in finite-difference methods, it is believed that this integral 
method, capable of producing solutions in reasonable agreement with a good 
finite-difference code in far less computer time, is a powerful tool for 
engineering calculations. Predictions of the integral method also agree well 
with experiment. The advantage of greatly reduced computation time, along 
with considerable robustness, makes the integral boundary-layer method well 
suited for use in interaction schemes. 


INVISCID ANALYSES 

During the course of this study, the boundary- layer method described in 
the previous section was interacted with a number of inviscid analysis codes. 

The three codes primarily used were of the FLO series, authored mainly by 
Jameson and Caughey (refs. 18 to 21). All three solve the full-potential 
equation, but in different forms or in different coordinate systems. This 
study was initiated using FLO-22, written for "wing only" configurations; no 
results for interactions using FLO-22 are presented in this paper. FLO-27, a 
wing-cylinder code, and FLO-30, a more general wing-body analysis code, were 
also used. Most of the interacted results presented here have been obtained 
using FLO-30. 

FLO-22 is the oldest code and, thus, is somewhat more well-tried than the 
others. The boundary-conforming grid used in FLO-22 is produced by shearing out 
the wing sweep, then applying a parabolic mapping, and finally shearing out the 
wing thickness. A planar vortex sheet is assumed. The full- potential equation 
is written in a nonconservative finite-difference form, using Jameson's "locally 
rotated” (ref. 18) operators. Solution is by line successive over-relaxation 
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(ISOR) with user-preset grid halving. As with all potential equation methods, 
good approximation is expected only in flows with weak shocks. Also, the non- 
conservative formulation produces solutions which do not maintain mass conserva- 
tion through a shock. The effect of the nonconservative form is. usually seen in 
the incorrectly predicted shock position. 

FLO-27 uses a finite-volume scheme to produce conservative differencing of 
the full-potential equation. After a preliminary Joukowsky transformation in 
the cross-flcw plane, which maps the cylindrical fuselage onto a vertical slit, 
the grid is generated in a manner similar to that in PLO-22, although another 
shearing transformation is applied to remove wing taper to maintain a constant 
number of grid points at each span station. Solution is again by LSOR, with 
additional iterative damping in the form of the mixed space-time derivative 
added to stabilize the scheme. Although the difference form used in FLO-27 is 
a better approximation than that in FLO-22, convergence is somewhat slower. 

Also, the finest grid used in FLO-27 is coarser than that in FLO-22 although 
FLO-27 requires more core storage. In all, FLO-27 was found to be about 1 0 per- 
cent more expensive to run than PLO-22 in these interaction studies. 

Most of the interaction logic was developed and initial operating experi- 
ence was gained by using FLO-27 for inviscid analysis. As shown later, however, 
the effect of the body on the inviscid flow is large for practical configura- 
tions. This necessitated the application of the developed interaction logic to 
a wing-body analysis code; FLO-30 was obtained for this purpose. The differenc- 
ing and solution of the full-potential equation as used in FLO-30 is quite simi- 
lar to that in FLO-27, in which a finite-volume scheme in conservative form is 
used. The coordinate system used, described in reference 19, is termed the 
"wind-tunnel" grid. Coordinates in the cross-flow plane are cylindrical, normal- 
ized with the fuselage section variation sheared out. A conformal mapping is 
used to "unwrap" constant-radius surfaces about the wing/wake surface similar 
to parabolic coordinates. As a last step in the grid-generation procedure, 
the wing thickness is sheared out. As with all finite-volume coordinate sys- 
tems, grid extent is finite. 

In order to model the effect of the boundary layer on the outer inviscid 
flow, the displacement thickness computed by the boundary- layer method was added 
to the basic wing shape in the surface normal direction. This technique was 
used, as opposed to using surface transpiration boundary conditions, for simplic- 
ity of incorporation into the FLO codes. The techniques are mathematically 
equivalent to (refs. 22 and 23), although some researchers have 

reported differences in results in practice (ref. 5) . The displac^nent thick- 
ness distribution is under-relaxed between iterations to prevent oscillations in 
total lift and lift distributions in the inviscid calculations. In addition, 
chordwise constant- value extrapolation is sometimes used on the first iteration 
when separation is predicted in the cove region of the lower surface; this was 
found to speed convergence greatly in difficult cases. 

The wake in the FLO codes is modeled as a contact discontinuity, across 
which the pressure is continuous but tangential velocity discontinuous. The 
small-disturbance form of this condition leads to a constant potential jump 
being carried downstream from the trailing edge on constant-span lines, with 
the value of the potential jump determined by using potential values at the 
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trailing edge. An approximate wake center-line shape is used, dependent on the 
wing trailing-edge angle, and is sheared out along with the wing surface shape 
during the coordinate generation process. This procedure fixes the wake along 
a coordinate surface to simplify the application of the wake boundary conditions. 


WAKE TREATMENT 

One major purpose of this work was to study the effect of the wake approxi- 
mation used. A full, asymptotically correct wake treatment for three-dimensional 
wing configurations of the kind described in reference 1 for two dimensions is 
at the present time lacking. However, displacement and curvature were the first 
wake effects accounted for in two-dimensional studies (ref. 22); therefore, a 
similar treatment is given here. 

For this study, the wake model used in the original FLO-30 code was 
replaced with one satisfying flew tangency on the wake displacement body and the 
pressure jump condition from wake curvature. As on the wing surface, the solid 
displacement surface model was used for the wake, constructed in two-dimensional 
strips. The pressure jump across the wake displacement body, related to the 
streamline curvature and momentum thickness of the wake, was reformed into a 
condition on the change in potential jump along the wake, as in reference 1. 

A fixed approximation to the wake center-line location was used for sim- 
plicity. This surface leaves the trailing edge smoothly at the averaged local 
trailing-edge angle, and the angle between the wake center-line surface and free- 
stream decays logarithmically as the general streamline shape of a point vortex 
in uniform flow. Wake displacement thickness was then applied along this 
surface. 

The flow curvature at the actual wake was computed from the rate of change 
of the. flow angle along the fixed wake; that is. 


de 

ds 


w 


Here, as in reference 1 , the approximation is made that the actual wake lies 
close to the preset inviscid wake near the trailing edge where the curvature 
effect is significant. An alternate approximation from inviscid theory was 
tried: 


h du^ 


(u dz> 


w 


where z is the direction normal to the wake. However, the cell spacing in the 
normal direction was found to be too large, and the curvature computed from the 


expression was too lew near the trailing edge. This being the region where wake 
curvature effects are largest/ the normal derivative expression for the curva- 
ture was dropped. Following reference 1 , the wake potential jump condition is 


dTw 
ds„ “ 


with = E'^w= ^u " % wake (where the subscripts u and S, indi- 

cate upper and lower, respectively) . 

The integral properties of the wake were computed during the same stage of 
the interaction Iterations as the boundary-layer calculations. A conpresslble 
lag-entrainment method, outlined in reference 13, was used to compute and 

0^, on streamwlse strips, with boundary-layer conditions along the trailing edge 
used as initial conditions. Two separate wake computations were made at each 
span station, for either side of the wake center line; the total integral prop- 
erties of the wake are then the sum of the corresponding properties from both 
sides. 


INTERACTION PROCEDURE 

Viscous-inviscid interaction results were produced in this study by iterat- 
ing potential, boundary-layer, and wake solutions until overall convergence was 
obtained. A cycle, or "global iteration" of this process, consists of comput- 
ing the external inviscid flow over the displacement surface of the wing and 
wake, with the inviscid wake potential jump condition modified by the curvature 
condition described previously, and, subsequently, of updating the viscous param- 
eters by boundary-layer and wake calculations, using the external velocity dis- 
tributions previously computed. The updated viscous parameters then provide a 
new displacement surface for the inviscid calculation Of the next global itera- 
tion. In practice, the displacement surface updates must be under relaxed with 
the viscous parameters of prior cycles to minimize oscillations. The computed 
displacement thickness on the wing is smoothed to give the displacement surface 
continuous curvature. Also, the inviscid calculations within each global itera- 
tion are stopped short of convergence within the potential equation iterative 
solution scheme. This practice was found to reduce the overall work involved 
in achieving global convergence. 

The potential codes described previously require very large amounts of 
computer storage. The finest grids used in FLO-22, FLO-27, and FLO-30 were 
(192 X 24 X 32), (1 61 x 18 x 35), and (1 61 x 24 x 32), respectively. Because 
of this, the integral boundary- layer code and interpolating codes required to 
process velocity and displacenent thickness data were kept as separate routines. 
The interaction iterations and flow of information between these programs were 
controlled by a job stream using the powerful CDC CYBER Control Language (CCL) . 
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RESULTS 


Experimental Configurations 

In order to show the advantages of the present more complete wake treat- 
mentf as opposed to the treatments previously used for three-dimensional wing 
calculations, configurations for which viscous effects are large must be used. 
The first test case chosen was selected from a family of supercritical wings 
tested at the Langley Research Center for advanced transport configurations 
(ref. 24) . The specific wing-body configuration used here is also one of the 
cases given in reference 4. A photograph of the wind-tunnel model, mounted for 
testing in the Langley 8-Foot Transonic Pressure Tunnel, is shown in figure 6, 
and a line drawing of the configuration is given in figure 7(a). The wing is 
of relatively high-aspect ratio and is swept and tapered. Supercritical air- 
foil sections with large aft camber are used over the span; a small amount of 
twist is employed, and the thickness ratio varies from 0.144 at the root to 
0.106 at the tip. The wing is mounted in a low position on a relatively wide 
body. 


In the second test case, a somewhat lower aspect-ratio wing is mid-mounted 
on a cylindrical fuselage. A schonatic of this configuration is shown in fig- 
ure 7(b) . Viscous effects for this case are less than on the first test case 
configuration, but a shock of greater strength appears on the upper surface for 
simil6u: free-stream Mach number eind overall lift level. The same airfoil sec- 
tion is used throughout the span, and the wing is slightly twisted. Data for 
this case were presented in references 25 and 26; this configuration is denoted 
wing A mid-mounted in those references. 


Basic Ccmputational Results 

The flow about the first configuration was computed by using FLO-30 alone 
(i.e., inviscid calculation) and using FLO-30 interacted with the integral 
boundary-layer scheme, with and without the viscous wake model. The cases were 
run at the experimental angle of attack of 1.65*^ and Mach number of 0.78; the 
Reynolds number was 2.4 x 10^ based on the meeui chord. 

Computed chordwise pressure distributions for the first test case configura- 
tion, taken at 24, 64, and 82 percent semispan, are compared with experiment in 
figure 8. The experimental stations are 26, 63, and 82 percent semispan. Agree- 
ment again is good, with cove-region and trailing-edge pressure predictions 
quite reasonable. Discrepancies at the upper surface leading edges of the out- 
board stations are noted; no explanation can be given at this time. 

Isoinclines of 6* for the first test case are shown in figure 9 for the 
upper and lower surfaces, respectively. Note that on the upper surface, lines 
of constant 6* are for the most part parallel to the leading edge, indicating 
little cross flow in the boundary layer. Near midspan, however, 6 grows more 
rapidly as the trailing edge is approached than in the inboard region. This 
strongly reduces the trailing-edge angle of the displacement body, effectively 
decamber ing this region more than inboard. Isoinclines of 6* on the lower 
surface are seen to sweep forward from parallel to the leading edge after mid- 
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chord, indicating more rapid relative thickening of the boundary layer toward 
the tip. In fact, the maximum 6*/c in the lower surface cove region is more 
than 30 percent greater in the outboard region than inboard. This reflects the 
cross flow in the boundary layer causing mass to build up toward the tip; the 
outboard sections are thus decamber ed more by viscous effects. 

In figure 1 0 is shown a ccnparison between the spanwise section lift dis- 
tributions from a purely inviscid calculation euid the interacted calculation 
using the curved wake model outlined in the previous section. FLO-30 was the 
inviscid method used in both calculations. Also shown in the figure are experi- 
mental values of c^ obtained from integrated pressure distributions. Apparent 
from the figure is the large decrease in lift due to the decambering effect of 
the boundary layer. This decambering is stronger toward the tip than near the 
root. A large amount of spanwise flow occurred in the boundary layer in the 
lower surface cove region; therefore, this effect was expected. The interaction 
method with the full-wake model agrees well with experiment, with a small dis- 
crepancy near the body. This difference is probably due to incorrect modeling 
of the flow in the wing-body juncture region. The total lift coefficient com- 
puted for the wing-body configuration is also very close to that from experimen- 
tal balance data - 0.472 computed versus 0.471 measured. 

Flow about the second test case configuration was also computed by using 
the interaction procedure with FLO-30, the integral boundary-layer scheme, and 
the full wake model. Again, experimental free-stream conditions were used; 
the Mach number was 0.819, the angle of attack was 1.96°, and the Reynolds num- 
ber was 6 >« 1 0® based on the mean chord. Transition was fixed at 5 percent 
chord. 


In figure 11 is shown a ccmparlson between chordwise pressure distribu- 
tions from the experimental data of references 25 and 26 and that computed with 
the interaction procedure. Shock strength and position are predicted well over 
most of the span, with a disagreement in position of a few percent chord develop- 
ing toward the tip. Cove region pressures on the lower surface are also pre- 
dicted reasonably well. Predicted trailing-edge pressures are slightly too high. 
However, the airfoil sections used in this wing have sharp trailing edges; 
therefore, this discrepancy cannot be attributed to inadequate modeling of blunt 
trailing-edge flow. These results are perhaps a good case for the development 
of an adequate three-dimensional strong interaction theory for wing trailing 
edges. 

Viscous effects for the wing A test case are also fairly large, as can be 
seen in the comparison of interacted and inviscid chordwise pressure distribu- 
tions at 68 percent semispan shown in figure 1.2. Without viscous correction, 
FLO-30 predicts the shock 0.1 4c too far aft cuid predicts a much greater cove 
region pressure. 


Wake Effects 

The results of calculations including wake effects were compared with those 
using the wake model of reference 19. The wake model in that procedure is modi- 
fied from that Incorporated in the original FLO codes only in the carrying of 



the trailing-edge thickness downstream to the edge of the computational domain 
along the assumed wake center line* Flow tangency is not strictly satisfied at 
the wakef and the pressures on either side of the wake are set equal only to the 
order of small-disturbance theory. 

In figure 13/ the parameter R is plotted against span distance, where 


/inviscid “ ^Z 

R E 

^Z f inviscid 


Thus, R is the normalized lift reduction due to viscous effects. The two 
curves shown in figure 13(a) are computed by using the present full wake model 
and by using the model of reference 19. Also shown are experimental values of 
R, computed using the experimental c^^s and the computed inviscid c^'s. Wake 
effects cire thus seen to be very important in the calculation of the lift dis- 
tribution on three-dimensional wings. In figure 13(b) the spanwise distribu- 
tions of the lift reduction parameter R computed using the present wake model 
are compared with those computed using the displacement thickness of the wake 
only, that is, neglecting the wake curvature effects. The experimental points 
are again shown. Note that the wake curvature effects are strongest again near 
the tip and that agreement in this area is much improved by including these 
effects. The curvature effect of the wake is seen to be larger than the dis- 
placement effect, which is similar to the results shown in two dimensions in 
reference 1 . 

In figure 14 are compared the chordwlse pressure distributions near mid- 
span from calculations using the full and original wake models. Note that the 
cove region pressure predicted using the full wake model is lower than that of 
the original model, in better agreement with experiment. 


Body Effects 

An early, unpublished version of the present interaction procedure used 
the full- potential code FLO-27 in wing-alone mode and the original wake model 
coupled with the integral boundary-layer method. Gross discrepancies between 
predicted and experimental pressure and lift distributions were seen by using 
this method. A comparison between lift distributions predicted using FLO-30 
with the original wake model and this early method for the first test case con- 
figuration is given in figure 15. The early method was run with a Mach number 
shift of 0.005 in an attempt to account for the body thickness effect. The 
test wing, low-mounted on a wide body, lies in a strong downwash field, which 
affects the predicted lift distribution significantly; less than a degree of 
downwash would be required to produce the lift discrepancy seen in figure 15. 
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Canputational Efficiency 


In all interaction procedures developed in this study/ it was found that 
the additional computational work required to converge a viscous-inviscid 
interaction case over that needed to converge the same case using the inviscid 
procedure alone was comparatively small. Also, inclusion of the wake effects 
was found to give no convergence penalty in the interaction procedure. With 
the full wake model, in fact, by properly selecting the S* under relaxation 
factor and the scheduling of the number of inviscid fine-grid iterations between 
viscous updates, the total number of fine-grid iterations required to converge 
an interacted case was only slightly greater than that required to converge a 
purely inviscid case to the same residual level. To converge to a level where 
the maximum residual in the field was 1 .5 x 10“®, the purely inviscid run 
required the work equivalent of 728 fine-grid Iterations on a halving sequence 
of three grids. The interacted test case required the equivalent of 756 fine- 
grid iterations. For the first two inviscid computations in the interaction 
sequence, calculation was begun afresh on the coarsest grid. Starting with the 
third iteration, the solution was restarted on the finest grid and a relatively 
small number of fine-grid iterations were run per viscous update. A total of 
seven global iterations were required, and a 6* underrelaxation factor of 0.6 
was used. Another interesting statistic is that, of the total CP time spent 
converging an interacted case, only 3 percent is spent computing the viscous 
parameters by using the integral method, compared with over 50 percent for the 
method of reference 6. 


CONCLUDING REMARKS 

A viscous-inviscid interaction computational method has been developed for 
three-dimensional transonic wing-body configurations. The procedure employs a 
three-dimensional integral boundary-layer method for viscous correction on the 
wing surfaces which produces results in good agreement with a finite-difference 
method in a fraction of the computer time. The integral method is stable and 
robust and incorporates a model for computation in a small region of streamwlse 
separation. A locally two-dimensional wake model, accounting for thickness and 
curvature effects, is also included in the interaction procedure. Computation 
time spent in converging an interacted result is, many times, only a small 
amount greater than that required to converge an inviscid calculation to the 
same maximum residual. 

Results for test cases in which viscous effects are large have shown good 
agreement with experimental data. In particular, the use of the full-wake model 
improved the prediction of spanwise load distribution and cove region pressure 
over results by using the wake model originally employed in PLO-30. For tran- 
sonic wing-body configurations, it is also apparent that body effects may be 
large. It is seen that a three-dimensional boundary- layer method must be used 
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oh the wing, rather than a two-dimensional strip method, in order to properly 
predict the increased decambering of the sections near the tip, due to cross 
flow in the boundary layer. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 29, 1981 
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Figure 1.- Breakdown of three- 
dimensional boundary layer 
into components. 
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Figure 3.- Comparison of displacement Figure 4.- Comparison of displacement 
thickness predictions for f) = 0.6. thickness predictions for T| s 0.7. 
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Figure 5.- Comparison with experiment for 
two-dimensional boundary layer only. 
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Figure 12.- Ccnparison of chordwise pressure 
distributions for wing A test configuration, 
n = 0.68. 
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Figure 13.- Lift reduction parameter plotted against span for advanced 

transport configuration. 
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Figure 14.- Chordwise pressure distributions 
for advanced transport configuration, 
n = 0.56. 
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Figure 15.- Influence of fuselage on spanwise lift dis 
tributions for advanced transport configuration. 
Boundary layer interaction only. 
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